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L  INTRODUCTION. 

This  report  presents  two  great  circle  navigation  algorithms,  a  closest  point  of  approach 
(CPA)  algorithm  and  an  intercept  algorithm.  It  also  presents  an  implementation  program 
that  is  written  in  the  BASIC  language  for  an  IBM  PC.  Instead  of  using  classical  spherical 
geometry  or  the  general  spherical  triangle  [Ref.  1],  these  algorithms  incorporate  rectangular 
coordinates  and  vectors  on  the  surface  of  the  sphere.  Portions  of  the  vector  formalism  were 
suggested  by  Reference  2. 

The  intent  of  the  report  is  to  provide  algorithms  for  spherical  earth  models  that  can 
be  used  for  situations  in  which  flat  earth  models  are  not  appropriate,  but  that  do  not 
require  the  sophistication  of  rotating  earth  models. 

n.  RECTANGULAR  COORDINATES  AND  VECTORS  ON  A  SPHERE 

In  a  spherical  earth  model,  a  point  P  on  the  earth's  surface  can  be  located  by  a  position 
vector  p  =  xi+ yj+zk  in  a  rectangular  coordinate  system  with  origin  at  the  earth's  center. 
In  matrix  form,  and  with  i,  y  and  z  expressed  in  spherical  coordinates, 


P  = 


r cos ^ cos  A 
rcos^sinA 

rsin</> 


where  <j>  is  the  latitude  and  A  is  the  longitude  at  the  point  and  r  is  the  distance  of  the 
point  from  the  earth's  center.  See  Figure  1. 
In  terms  of  the  unit  vector  x  where 


x  = 


XI 

cos  <f>  cos  A 

zj 

= 

cos  <f>  sin  A 

x* 

sin0 

(1) 


p  =  rx.  We  can  think  of  x  as  the  unit  vector  normal  to  the  surface  at  the  point.  It  is 
convenient  to  define  two  other  unit  vectors  in  the  tangent  plane  to  the  earth's  surface  at 
the  point.  These  are  local  north  n  and  local  east  e,  defined  by 


x<a 
n  = 


X4 


where 


sin  <f>  cos  A 
sin  <f>  sin  A 

C08(j> 


and 


e  = 


sin  A 
cos  A 
0 


(2) 


and 


*A  = 


fa 


The  vectors  x,  n  and  e  provide  the  basis  for  a  right-handed  orthogonal  coordinate  system. 
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Figure  1.  The  Earth  Coordinate  System 


Figure  2.  The  Course  Vector 


Let  a  course  a  be  designated  at  a  point  on  the  earth's  surface.  We  wish  to  determine 
the  unit  course  vector  c  in  the  direction  a  which  lies  in  the  tangent  plane  at  the  point 
in  terms  of  <p,  A  and  a  (see  Fig.  2).  An  arbitrary  vector  a,  defined  at  the  point  P  with 
coordinates  (r,  A,  <j>)  and  lying  in  the  tangent  plane  at  the  point,  can  be  rotated  clockwise 
(looking  from  P  toward  the  origin)  through  an  angle  a  by  using  the  operation  Rp(a)a, 
where  Rp(a)  is  the  rotation  operator.  Rp(a)  is  a  composite  rotation  in  3-space  and  can  be 
decomposed  into  fundamental  rotations  in  one  of  several  ways.  One  way  is  to  proceed  as 
follows:  First,  move  P  (carrying  a  with  P)  along  its  parallel  of  latitude  to  the  x-z  plane  (the 
Greenwich  meridian);  this  is  equivalent  to  a  clockwise  rotation  about  the  2- axis  through 


an  angle  A  and  is  denoted  by  R»(A).  Next,  move  P  along  the  Greenwich  meridian  to 
the  equator:  this  is  equivalent  to  a  counterclockwise  rotation  about  the  y-axis  through  an 
angle  <f>  and  is  denoted  by  Ry(-0).  The  point  P  now  lies  on  the  x-axis  with  coordinates 
(r,  0, 0)  and  the  vector  a  at  P  makes  the  same  angle  with  the  Greenwich  meridian  as  it  did 
with  respect  to  the  original  meridian  at  (r,  A,  <j>).  Next,  leave  P  on  the  Greenwich  meridian 
at  the  equator  and  rotate  a  through  an  angle  a  about  the  i-axis;  this  clockwise  rotation 
is  denoted  by  Rx(a).  The  vector  a  has  now  been  rotated  through  the  desired  angle  a 
with  respect  to  the  Greenwich  meridian.  When  P  is  returned  to  its  original  position  by 
reversing  the  steps  which  got  it  to  the  equator  on  the  Greenwich  meridian,  it  will  have 
been  rotated  through  the  angle  a  with  respect  to  the  original  meridian  of  Pt  i.e.  Ry(0) 
followed  by  R«(-A).  The  composite  rotation  operator  Rp(a)  is 

Rp(a)  =  Ra(-A)RyMRs(cOR,M)R,(A). 


The  course  vector  can  then  be  written  as  c  =  Rp(ot)n.  The  fundamental  z-,  y-  and  z-axis 
rotation  operators  are  given  by 


RxM  = 


Ry(9)  = 


R*M  = 


10  0 

0       cosl    sin0 
0     —  sin  9    cos  9 

cos0    0    -sin0 
0       1        0 

sin  9     0        cos  9 

cos#     sinfl     0 

-  sin  9    cos  9    0 

0  0       1 


and 


These  rotation  operators  are  consistent  with  a  right-handed  coordinate  system  and  positive 
signs  of  9  for  a  counterclockwise  rotation  of  the  coordinate  system  or  a  clockwise  rotation 
of  a  point  P  about  the  x-,  y-  or  z-axes,  respectively,  as  viewed  looking  toward  the  origin 
from  the  positive  end  of  the  rotation  axis  [Ref.  3,  pg.  43  and  Ref.  4,  pg.  100].  Some 
simplification  in  determining  Rp(a)n  can  be  obtained  by  noting  that 


R„  M)Rz(A)n  = 


=  k. 


Thus  c  can  be  found  from 


c  =  Rl(-A)R,(pI(a)k. 
3 


(3) 


If  an  object's  course  vector  c  is  known  at  some  point  with  position  vector  p  then  it  is 

easily  shown  (see  Fig.  3)  that 

n  •  c  =  cos  a    and 

e  •  c  =  sin  a. 

From  these  relations,  the  course  a  is  found  to  be 

a  =  qatn(e-c,n-e)  (4) 

where  qatn  is  a  quadrant  determining  arctangent  function  (see  Appendix  A  and  Ref.  5). 


e 


Figure  3.  An  Object's  Course 


The  course  can  also  be  determined  from  the  great  circle  vector  g  that  is  defined  by 

g  =  xxc,  (5) 

where  p  =  rx.  Prom  Figure  3  we  see  that  the  following  relationships  hold: 

n    g  =  cos  £  =  sin  a     and 
e  •  g  =  cos  r)  =  -  cos  a, 

whence 

a  =  qatn(n  g,  -e  •  g)  (6) 

If  the  object  is  traveling  with  speed  v  and  is  not  maneuvering,  it's  course  will  be  a 
great  circle.  Let  vo  =  vco  denote  the  object's  velocity  vector,  where  Co  is  its  course  and 
po  is  its  position  vector  at  time  0.  At  some  subsequent  time  t,  the  object's  position  vector 
will  be  p(i)  and 

p(t)  =  ap0  +  bv0t  (7) 


v0 

Po 

vt 

?t 

p 

vt/r     / 

Figure  4.  The  Velocity  Vector 

where  a  and  b  are  to  be  determined  (see  Fig.  4).  Dotting  Equ.  7  from  the  right  by  p0,  we 
see  that 

p(t)  •  po  =  ap0  •  po 

or,  with  angles  b  radians, 


_  P(0  •  Po  .     1    m   n   _  1 


o       fvf 

r cos  I  — 


or 


po  •  Po       r- 

fvt\ 
a  =  cos    —  ) . 

Similarly,  dotting  Equ.  7  from  the  right  by  vo  we  find 

p(0*vo  =  &vo  -To*=  btPt 


or 


so  that 


Thus, 


b  =  -«tp(0  •  vo  =  -or 


rvcos 


h  =  7t*m(j) 


p(i)=Pocos(M  +v0^sinfy 


In  terms  of  the  unit  vectors,  this  equation  becomes 

rx(t)  =  nco cos  I  ~  J  +  vc0^ sin  (  —  j  , 


or 

x(f )  =  xo  cos  I  —  j  +  c0  sin  (  —  J  .  (8) 

Applications  of  the  these  relations  are  made  in  the  following  sections  of  this  report. 

m.  GREAT  CIRCLE  NAVIGATION. 

The  aDirect  Solution  Algorithm*  computes  the  latitude  and  longitude  of  a  position  P2 
and  the  backward  azimuth  from  P2  to  Pi  given  the  latitude  and  longitude  of  a  position 
Pi,  the  forward  azimuth  from  Pi  to  Pj  and  the  distance  from  Pi  to  Pj.  The  "Inverse 
Solution  Algorithm"  computes  the  distance  from  position  Pi  to  position  Pg,  the  forward 
azimuth  from  Pi  to  P2,  and  the  backward  azimuth  from  P2  to  Pi  given  the  latitude  and 
longitude  of  positions  Pi  and  P2.  Details  of  these  algorithms  using  the  concept  of  the 
general  spherical  triangle  are  presented  in  Reference  5.  These  models  are  redeveloped  here 
using  the  concepts  of  Section  II. 

A.  The  Direct  Solution  Algorithm.  Given  Pi(^i,  Ai),  forward  azimuth  (bearing)  a^ 
and  distance  d,  find  <fo  and  A?  of  P2  and  the  backward  azimuth  c*2i-  Proceed  as  follows: 

1.  Prom  <f>\  and  Ai,  compute  the  components  of  Xi  using  Equ.  1. 

2.  Compute  Ci  from 

d  =  Ra(-Ai)IM0i)lUai2)k. 

3.  Compute  x2  using 

fd\  •    fd\ 

xj  =  xi  cos  1  -  I  +  ci  sin  1  -  ]  . 

Note,  with  d  =  vt  in  nautical  miles  and  r  =  60(180°/ir),  Equ.  8  becomes 


Xj  =  Xi  cos 


(m)+Ci™(£>) 


where  the  arguments  of  the  cosine  and  sine  are  now  in  degrees. 

4.  Prom  the  components  of  X2  compute 

<f>2  =  sin-1  (^2)     and 
A2  =  qatn(zj2,a:i2). 

5.  Compute  g  =  Xi  x  ci. 

6.  Compute  n2  and  e2  using  Equ.  2. 

7.  Using  Equ.  6  compute  0:21  =  -  qatn(n2  •  g,  -e2  •  g).  Note  that  the  sign  change  occurs 
because  a2i  is  the  backward  azimuth. 
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B.  The  Inverse  Solution  Algorithm.  Given  Pi  (fa ,  Ai )  and  P2(^2,  A2),  find  the  distance 
d  from  Pi  to  P2,  the  forward  azimuth  a  12  from  Pi  to  P2  and  the  backward  azimuth  ar2i 
from  P2  to  Pi .  Proceed  as  follows: 

1.  From  fa  and  At-  compute  x«,  for  i  =  1, 2. 

2.  Compute  d  =  rcos_1(xi  •  x2),  whence  d  =  60(180/ir) cos-1(xi  •  x2)  is  the  distance  in 
nautical  miles. 

3.  Compute  n«  and  e,-  for  t  =  1, 2  (see  Equ.  2). 

4.  Compute 

Xi  XX3 

g  = 


|xi  xx2| 

This  equation  arises  because  the  great  circle  passes  through  both  Xi  and  x2,  hence  g 
must  be  perpendicular  to  Xi  and  to  x2. 

5.  Compute  ari2  =  qatn(ni  •  g,  -ei  •  g). 

6.  Compute  or2i  =  -  qatn(n2  •  g,  -e2  •  g). 

IV.  CLOSEST  POINT  OF  APPROACH  (CPA)  PROBLEM. 

Consider  two  objects  traveling  on  different  great  circle  paths.  From  Equ.  8,  their  tracks 
will  be  characterized  by  the  equations 

x,(f)  =  x,o  cos  Uit  +  c,o  sin  w»f,    for  i  =  1,2,  (9) 

where  w,-  =  t/,/r.  At  any  time  t,  their  angular  separation  $(t)  in  radians  is  determined  by 

cosa(*)=xi(i)-x2(0.  (10) 

The  time  of  minimum  separation  (CPA)  is  that  time  T  when  ^[cosa(i)]  =  0.  That 
is,  we  must  find  T  such  that 


xl(r).^  +  ^.X3(r)  =  o. 


dt  dt 

Unfortunately,  this  equation  cannot  be  solved  explicitly.    One  approach  is  to  use  the 
Newton- Raphson  iteration  method  [Ref.  6]  to  find  the  root  T  of  /(*),  where 

f{t)  =  Xi (f  J  •  — ^—  +  — ^—  •  X2(*J. 

Taking  the  required  derivatives  of  Equ.  9  and  performing  the  required  dot  products,  we 
find  that 

f(t)  =  j4sinu;i*sinii;2*  +  #coswi*  cosw2<  +  Csinu/ifcoso;2f  +  £coswitsinw2*, 
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where 

A  =  — (wiXio  •  C20  +  W2X20  •  cio)i 

B  =  uj\ c10  •  X20  +  W2C20  '  xio» 

C  =  -(uiXio  •  X20  -  W2C20  •  cio)    and 

D  =  <jj\  C]q  •  C20  ~~  W2X20  •  Xio- 

The  use  of  the  Newton- Raphson  method  requires  that  the  derivative  of  /(f)  with  respect 
to  t,  namely  /'(f),  be  computed  or  estimated.  We  find  that 

f{t)  =  -  (CW2  +  .Dwi)sinwi£sina;2*  +  (Du^  +  Cui)  cos u\t cos u^t 
+  (Au>2  -  Bui)  sin. uit cos uzt  -  (BU2  -  Au^coswitsmurf. 

The  Newton- Raphson  method  requires  us  to  compute 

where  to  is  an  initial  estimate  of  T,  until  some  value  of  %  is  found  for  which  f(U)/f{U)  is 
sufficiently  close  to  zero. 

The  CPA  option  in  the  computer  program  (Appendix  B)  will  print  out  the  time  to 
CPA  (from  time  t  =  0),  the  distance  between  objects  at  CPA  and  the  bearing  from  object 
1  from  object  2  at  CPA.  Also  printed  is  the  number  of  iterations  required  for  convergence 
of  Equ.  11  to  \f(U)/f'{U)\  less  than  0.00001  hours.  A  negative  time  to  CPA  indicates  that 
CPA  has  already  occurred. 

V.  INTERCEPT  PROBLEM. 

The  intercept  problem  is  divided  into  two  problems,  each  of  which  may  require  an  answer. 
In  both  problems  we  are  given  the  initial  position,  course  and  speed  of  a  target  as  well 
as  the  position  of  an  interceptor  or  launch  platform.  In  the  first  problem  we  are  given 
the  time  (or  elapsed  time)  at  which  intercept  is  desired  and  are  required  to  compute  the 
vehicle  speed  needed  for  an  intercept  to  take  place.  In  the  second  problem  we  are  given 
the  speed  of  an  intercept  vehicle  and  wish  to  compute  the  time  required  for  an  intercept 
to  occur  provided  an  intercept  can  be  made.  Provision  for  both  of  these  options  is  made 
in  the  program  presented  in  Appendix  B. 

A.  Speed  Required  to  Intercept.  Inputs  are  the  initial  latitude  and  longitude  of  an 
interceptor  and  a  target,  and  the  target  course  and  speed.  Also  input  is  the  time  of 
the  desired  intercept.   Outputs  are  the  speed  required  of  the  interceptor,  the  course  of 
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the  interceptor,  the  distance  traveled  to  intercept,  and  the  latitude  and  longitude  of  the 
intercept. 

From  the  inputs,  use  Equ.  1  to  compute  xio  and  x^o,  the  position  vectors  of  the 
interceptor  and  target,  respectively,  at  time  t  =  0.  Denote  the  time  required  to  intercept 
by  t.  Compute  the  target  course  vector,  C20,  using  Equ.  3  and  then  compute  Xa($),  the 
position  of  the  target  at  the  time  of  intercept,  using  Equ.  8.  The  remainder  of  the  problem 
is  solved  using  the  inverse  solution  algorithm  discussed  previously.  The  speed  required  for 
intercept  is  given  by  t/i  =  d/t,  where  d  is  the  distance  from  the  initial  interceptor  position 
to  the  target  position  at  the  time  of  intercept. 

B.  Time  Required  to  Intercept.  Inputs  are  the  initial  latitude  and  longitude  of  an 
interceptor  and  a  target,  and  the  target  course  and  speed.  For  a  given  interceptor  speed, 
it  may  or  may  not  be  possible  to  make  an  intercept.  We  develop  two  sub- algorithms.  The 
first  algorithm  is  to  compute  the  minimum  interceptor  speed  required  to  achieve  intercept 
and  the  time  required  to  make  such  an  intercept.  The  second  algorithm  queries  the  user  to 
input  an  interceptor  speed.  K  the  speed  is  at  least  that  required  for  intercept,  then  the  time 
required  to  intercept  is  computed.  (If  the  interceptor  speed  is  greater  than  the  minimum 
required,  there  are  two  possible  solutions  for  the  time  to  intercept — the  shortest  time  to 
intercept  is  computed  by  the  algorithm).  Outputs  are  the  minimum  required  interceptor 
speed,  the  time  to  intercept  at  minimum  speed,  the  course  of  the  interceptor,  the  distance 
traveled  to  intercept,  and  the  latitude  and  longitude  of  the  intercept. 

The  first  problem  is  to  determine  the  minimum  speed,  vm,  required  to  make  an  in- 
terception. This  can  be  accomplished  by  finding  the  time  of  intercept  tm  which  makes 
dv/dt  =  0.  We  can  relate  v  to  a,  the  angular  separation  in  radians,  between  the  two  points 
xio  and  X2{t)  by  the  relation  v(t)  =  rs(t)/t.  We  find  that  dv(t)/dt  =  0  implies 


Ids      s 


=  0.  (12) 


Anticipating  that  it  will  not  be  possible  to  find  a  closed  form  solution,  we  prepare  to  use  the 
Newton- Raphson  procedure  (Equ.  11).  Multiplying  Equ.  (12)  by  t2  gives  us  the  function 

for  which  we  wish  to  find  tm  such  that  f{tm)  =  0.  Also  needed  is 

/'«)=<§• 
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Using  Equ.  10  to  determine  s(t)  we  find  that 


ds 

dt 

d*3_ 
dt* 


1 


sin  3 


sins 


-xio 


dx2{t) 


dt 


cos  a 
sin2s 


(xio 


and 


dx7(t)\ 


dt 


) 


+  xio 


dt* 


where 


xio  •      ?,      =  -^2  [*io  •  x2o  sin  w2*  -  Xio  •  c2o  cos  w2i]     and 


'  — 5f3 —  ""       2  ^ 


Xqq  co8u2t  +  Xio  -  c2o  sin  w2i] . 


*m- 


The  Newton-  Raphson  procedure  continues  until  tm  is  found,  then  vm  =  sr/tf 

The  second  problem  is  to  find  the  time  of  interception  ti  when  the  interceptor's  speed 
V\  is  given.  Once  more  the  Newton-Raphson  procedure  is  used.  As  before,  we  can  relate 
Vi  to  a(t),  the  angular  separation  in  radians,  between  two  points  Xio  and  x2(0  by  the 
relation  v(t)  =  rs(t)/t,  which  tells  us  that  we  must  require  wi  =  a(f  )/t.  That  is,  we  wish 
to  find  (/  for  which  s(ti)/t[  equals  the  constant  wi,  or  equivalently,  we  wish  to  find  ti  such 
that  f{ti)  -  0  where 


Also  needed  is 


/(0  =  j-«i. 


™-\{i-iy 


The  equation  for  da/dt  is  the  same  as  that  given  in  the  previous  paragraph.  The  remaining 
output  is  found  using  the  inverse  solution  algorithm  for  the  points  xi0  and  x2(i/). 
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VL  SAMPLE  PROBLEMS 

Master  Menu.  The  master  menu  for  algorithm  demonstration  program  is  shown  below. 

ALGORITHM  DEMO 

1)  DIRECT  SOLUTION 

2)  INVERSE  SOLUTION 

3)  FIND  CPA 

4)  SPEED  NEEDED  TO  INTERCEPT 

5)  TIME  NEEDED  TO  INTERCEPT 

6)  QUIT 

Problem  1.  Suppose  you  are  at  San  Francisco  (latitude  37° 47'  north  and  longitude 
122° 25'  west),  that  your  initial  course  is  260°  and  that  you  travel  a  distance  of  4000  n.  mi. 
What  is  your  final  position?  Select  Option  1  from  the  master  menu. 

DIRECT  SOLUTION 

1st  LATITUDE   DD.MMSS  (-S)  ?  37.47 

1st  LONGITUDE  DDD.MMSS  (-E)  ?  122.25 

INITIAL  COURSE  DDD.MMSS  ?  260 

DISTANCE  (n.  mi.)  ?  4000 

SPHERICAL  EARTH  DIRECT  SOLUTION 
2nd  LATITUDE     6°41.9' 
2nd  LONGITUDE  -172°00.7' 
BACK  AZIMUTH    51°35.9' 

PRESS  ANY  KEY  TO  CONTINUE 
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Problem  2.  Suppose  you  are  at  San  Francisco  (latitude  37°  47'  north  and  longitude 
122° 25'  west)  and  that  your  destination  is  Sydney,  Australia  (latitude  33° 51'  south  and 
longitude  151°  13'  east).  How  far  do  you  travel,  what  is  your  initial  course,  and  what  is  the 
backward  azimuth  from  Sydney  to  San  Francisco?  Select  Option  2  from  the  master  menu. 

INVERSE  SOLUTION 


1st  LATITUDE   DD.MMSS 

(-S) 

? 

37.47 

1st  LONGITUDE  DDD.MMSS 

(-E) 

7 

122.25 

2nd  LATITUDE   DD.MMSS 

(-S) 

? 

-33.51 

2nd  LONGITUDE  DDD.MMSS 

(-E) 

? 

-151.13 

SPHERICAL  INVERSE  SOLUTION 

DISTANCE  6446.3  n.mi. 

FORWARD  COURSE        240°18.9; 
BACK  COURSE  55°45.9' 

PRESS  ANY  KEY  TO  CONTINUE 
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Problem  3.  Suppose  an  observer  is  at  20°  north,  60°  west  traveling  on  a  course  of 
010°  at  a  speed  of  15  knots.  A  target  is  reported  to  be  at  34°  north,  50°  west  on  a  course 
of  220°  at  a  speed  of  300  knots.  Assuming  that  neither  observer  or  target  changes  course 
or  speed,  how  much  time  will  elapse  until  CPA  and  at  CPA  where  will  the  target  be  with 
respect  to  the  observer?  Select  Option  3  from  the  master  menu. 


FIND  CPA 

1st  LATITUDE   DD.MMSS  (-S) 

?  20 

1st  LONGITUDE  DDD.MMSS  (-E) 

?  60 

INITIAL  COURSE  DDD.MMSS 

?  10 

SPEED  (knots) 

?  15 

2nd  LATITUDE   DD.MMSS  (-S) 

?  34 

2nd  LONGITUDE  DDD.MMSS  (-E) 

?  50 

INITIAL  COURSE  DDD.MMSS 

?  220 

SPEED  (knots) 

?  300 

TIME  TO  CPA  -  3h09m48s 

DISTANCE  AT  CPA  =  67.03  n.mi. 
BEARING  AT  CPA    -  304°06.3/ 
NO.   ITERATIONS    =  3 

PRESS  ANY  KEY  TO  CONTINUE 
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As  an  additional  output,  a  table  of  observer  positions  and  target  positions  is  given  at 
CPA  and  six  equally  spaced  times  before  and  after  CPA. 


FIND  CPA 

TIME 

LAT  1 

LONG  1 

DISTANCE 

BEARING(l->2) 

00s 

20°  00.0' 

60°  00.0' 

994.34 

30°  18. 8' 

31m38s 

20°  07. 8' 

59°  58. 5' 

829.45 

29°31.6' 

Ih03ml6s 

20°  15. 6' 

59°  57.1' 

664.78 

28°  21. 6' 

lh34m£4s 

20°23.4' 

59°  55. 6' 

500.56 

26°26.3/ 

2h06m32s 

20°31.2' 

59°  54.1' 

337.43 

22°  39. 8' 

2h38ml0s 

20°  38. 9' 

59°  52. 7' 

178.42 

12°02.7' 

3h09m48s 

20°46.7' 

59°  51. 2' 

67.03 

304°  06. 3' 

3h41m27s 

20°  54. 5' 

59°49.7' 

178.43 

236°  10.0' 

4hl3m05s 

21°02.3' 

59°48.2' 

337.43 

225°  33. 2' 

4h44m43s 

21°  10.1' 

59°46.7' 

500.58 

221°47.3' 

5hl6m21s 

21°  17. 9' 

59°45.3' 

664.81 

219°52.7/ 

5h47m59s 

21°25.7' 

59°  43. 8' 

829.5 

218° 43. V 

6hl9m37s 

21°33.4/ 

59°42.3' 

994.41 

217°57.6/ 

PRESS  ANY  KEY  TO  CONTINUE 
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Problem  4.  Suppose  an  observer  at  20°  north,  60°  west  wishes  to  launch  an  inter- 
ceptor at  a  target  reported  to  be  at  34°  north,  50°  west  on  a  course  of  220°  at  a  speed  of 
600  knots.  If  interception  is  required  to  take  place  45  minutes  (2700  seconds)  after  launch, 
how  fast  must  the  interceptor  travel,  and  where  will  the  intercept  take  place?  (Assume 
that  the  target  does  not  change  course  or  speed.)  Select  Option  4  from  the  master  menu. 

SPEED  NEEDED  TO  INTERCEPT  (l->2) 

1st  LATITUDE      DD.MMSS  (-S)  ?  20 

1st  LONGITUDE  DDD.MMSS  (-E)  ?  60 

2nd  LATITUDE      DD.MMSS  (-S)  ?  34 

2nd  LONGITUDE  DDD.MMSS  (-E)  ?  50 

2nd  COURSE        DDD.MMSS  ?  220 

2nd  SPEED  (knots)  ?  600 

TIME  TO  INTERCEPT  (SECONDS)   ?  2700 

SPEED  REQUIRED  -  730.1  knots 

BEARING  TO  INTERCEPT      =  26°06.9' 

RANGE  TO  INTERCEPT  -  547.5  n. mi. 

LATITUDE  OF  INTERCEPT    =  28°08.0' 

LONGITUDE  OF  INTERCEPT  -  55°27.6' 

1)  CHANGE  TIME  OF  INTERCEPT 

2)  NEW  PROBLEM 

3)  MASTER  MENU 
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Problem  5.  As  in  the  previous  problem,  suppose  an  observer  at  20°  north,  60°  west 
wishes  to  launch  an  interceptor  at  a  target  reported  to  be  at  34°  north,  50°  west  on  a 
course  of  220°  at  a  speed  of  600  knots.  If  the  maximum  speed  of  the  interceptor  is  700 
knots,  how  long  will  it  take  before  interception  can  occur,  what  should  be  the  interceptor's 
initial  great  circle  heading,  and  where  will  the  intercept  take  place?  (Assume  that  the 
target  does  not  change  course  or  speed.)  Select  Option  5  from  the  master  menu. 

TIME  NEEDED  TO  INTERCEPT  (l->2) 

1st  LATITUDE      DD.MMSS  (-S)  ?  20 

1st  LONGITUDE  DDD.MMSS  (-E)   ?  60 

2nd  LATITUDE      DD.MMSS  (-S)  ?  34 

2nd  LONGITUDE  DDD.MMSS  (-E)  ?  50 
2nd  COURSE        DDD.MMSS  ?  220 

2nd  SPEED  (knots)  ?  600 

MINIMUM  SPEED  REQUIRED  TO  INTERCEPT      =     52.6  knots 
TIME  REQ'D  TO  INTERCEPT  AT  MIN  SPEED    -     Ih39m50s 

INTERCEPTOR  SPEED  (knots)     ?  700 

TIME  REQUIRED  -  46m03s 

BEARING  TO  INTERCEPT  =  25°56.l' 

RANGE  TO  INTERCEPT  -  537.2  n.ml. 

LATITUDE  OF  INTERCEPT  =  27°59.6# 

LONGITUDE  OF  INTERCEPT  -  55° 34. V 

1)  CHANGE  INTERCEPTOR  SPEED 

2)  NEW  PROBLEM 

3)  MASTER  MENU 
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APPENDIX  A:  The  QATN  Function 

This  routine  is  the  standard  arctangent  function  corrected  for  quadrant.  The  quadrant  arc- 
tangent function  is  occasionally  implemented  as  the  ATAN2  function,  the  ANGLE  function 
or  the  Rectangular-to-Polar  function. 

Entering  variables  are  the  x-  and  y-coordinates,  X  and  Y.  The  exiting  variable  is 
the  angle  8,  where  —  ir  <  8  <  n.  Use  of  the  quadrant  arctangent  function  is  denoted  by 
8  =  qatn(y,X). 

1.  If  X  ■£  0,  go  to  step  4. 

2.  Set8  =  (*/2)*sgn(y). 

3.  Go  to  step  8. 

4.  Set  8  =  arctan(y/X). 

5.  If  X  >  0,  go  to  step  8. 

6.  Set  8  =  8  +  *  *  sgn(y). 

7.  If  y  =  0,  set  8  =  u\ 

8.  Return. 

Note: 

If  y  >  0  then  sgn(y)  =  +1. 
If  y  =  0  then  sgn(y)  =  0. 
Ify  <  0  then  sgn(y)  = -1. 

Users  of  Microsoft  BASIC  can  simplify  the  qatn  function  significantly  by  using  the  code 
given  below.  To  return  an  angle  of  8  (designated  by  A  in  the  code)  in  the  range  of  (-ff,  ir)} 
use: 

PI  »  4*ATN(1):     TP  =  PI  +  PI:  EPS  »  1E-33 

A  -  ATN(Y/(X-EPS*(X=0)))   -  PI*(X<0)*(SGN(Y)   -   (Y=0)) 

To  return  a  value  of  A  in  the  range  of  (0,  2tt),  use: 

PI  -  4*ATN(1):     TP  -  PI  +  PI:  EPS  »  1E-33 

A  =  ATN(Y/(X-EPS*(X=0)))   -  PI*(X<0)  +  TP*(X  >=  0)*(Y<0) 
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APPENDIX  B:  Program  Listing 

10  '  "RECT  COORD  ALGORITHMS-  R.H.  SHUDDE,  03-03-86.   REV.  03-19-86  1000 

13  '  RECTALGR-.A 

100  DEFDBL  A-Z 

110  PI4-ATN(l#):PI2-PI4+PI4:PI-PI2+PI2:TP-PI+PI:RD-PI/180#:EPS-lD-33 

120  AE=60#*360#/TP 

130  DEF  FNM(X)-X-MO*INT(X/MO) : '  X  MOD  MO  FUNCTION. 

140  DEF  FNL(X)-X-TP*INT((X+PI)/TP):'  LONGITUDE  ADJUST  (-PI.PI) 

150  DEF  FNR(X)-INT(X*M0+.5)/M0:'  ROUNDING  FUNCTION. 

160  DEF  FNG(X)=X+PI*SGN(X)*(AB8(X)>PI2):'  LATITUDE  ADJU8T  (-PI/2, PI/2) 

170  DEF  FNAC8(X)-ATN(8qR(l#-X*X)/(X-EP8*(XK)#)))-PI*(X<0#):,  ARCCOS 

180  DEF  FNASN(X)=ATN(X/(8QR(1#-X*X)-EPS*(ABS(X)=1#))):'  ARCSIN 

190  'QATN  (-PI.PI)  FUNCTION: 

200  DEF  FNATN2(Y,X)=ATN(Y/(X-EP8*(X=0#)))-PI*(X<0#)*(SGN(Y)-(Y=0#)) 

210  'QATN  (O.TWOPI)  FUNCTION: 

220  DEF  FNATNP(Y,X)=ATN(Y/(X-EP8*(X=O#)))-PI*(X<0#)+TP*(X>*0#)*(Y<0#) 

230  'CROSS  PRODUCT:  X( . )-Xl ( . )  X  X2( . ) : 

240  DEF  FNCX(X1,Y1,Z1,X2,Y2,Z2)-Y1*Z2-Z1*Y2 

260  DEF  FNCY(X1,Y1,Z1,X2IY2IZ2)-X2*Z1-Z2*X1 

260  DEF  FNCZ(X1,Y1,Z1,X2,Y2,Z2)»X1*Y2-Y1*X2 

270  GOTO  2000 

280  ' 

1000  '  DECIMAL  TO  HH  MM  SS 

1010  V$--  ":IF  X<0  THEN  V$-»-":X— X 

1020  X=X+l/720O#:Y=INT(X) :Z=LEN(8TR$(Y))-1 

1030  KX-0:IF  Y<>0  THEN  V$-V$+RIGHT*("  "+STR$(Y),Z)+"h":KX-l 

1040  X=60#*(X-Y):Y=INT(X) 

1050  IF  Y<>0  OR  KX-1  THEN  X|-STR$(100#*Y)  :V$-V$♦RIGHT$(X$,2)♦■n,, 

1060  X=60f  *(X-Y) :Y=INT(X) :X|=8TR$(100#+Y) : V$-V$+RIGHT$(X$, 2) +"a" : RETURN 

1070  ' 

1080  '  DECIMAL  TO  DDD  MM.F 

1090  V$-"  ":IF  X<0  THEN  V$-"-":X— X 

1100  X*X+1/1200#:Y-INT(X):V$-V$+RIGHT$("  "+STR$(Y),3)+CHR$(248) 

1110  X-600#*(X-Y) :Y-INT(X) :X$-8TR$(1000+Y) 

1120  V$=V$+MID$(X$.3,2)+V"+RIGHT$(X$,1)+B*":RETURN 

1130  ■ 

1140  ■  DDD.MMSS  TO  DECIMAL 

1150  IX-0:FOR  Zl-l  TO  LEN(V$) :C$-MID$(V$,Z! ,1) :IF  C$-"."THEN  IX-ZI 

1160  NEXT:  IF  IX«0  THEN  X«VAL(V|) :  RETURN 

1170  X-VAL(LEFT$(V$,IX)):SN-1:IF  X<0#  THEN  SN— 8N:X— X 

1180  V|»V$+"0000":Y=VAL(MID$(V$fn+l,2)):Z»VAL(MID$(V$,IX+3,2)) 

1 190  X-SN* ( (Z/60#*Y) /60#+X) :  RETURN 

1200  ' 

1300  'DIRECT  SOLUTION,  SPHER  EARTH-RECT  COORD.  ALL  ANGLES  MUST  BE  IN  RADIANS 

1310  'INPUT:  LATITUDE  PI.  LONGITUDE  LI.  FORWARD  AZIMUTH  A12  AND 

1320  '  DISTANCE  DD  TO  A  POINT  P2.  NOTE:  DD  IS  IN  RADIANS. 

1330  'OUTPUT:  LATITUDE  P2.  LONGITUDE  L2  AND  BACKWARD  AZIMUTH  A21. 

1340  P=P1:L=-L1:G0SUB  1750: 'CHNG  SIGN  OF  LI  GIVES  RIGHT-HANDED  COORDS 

1360  FOR  I!=l  TO  3:P081(I!)=P0SI(I!) :N0RTH1(I!)=N0RTH(I!) :EAST1(I! )=EAST(I!) : 

NEXT 

1360  GOSUB  1810:CD=COS(DD):SD=SIN(DD) 
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1370 
1380 
1390 
1400 
1410 
1420 
1500 
1510 
1620 
1530 
1540 
1550 

1560 
1570 

1580 
1590 
1600 
1610 
1620 

1630 

1640 
1660 
1700 
1710 
1720 
1730 
1740 
1750 
1760 
1770 
1780 
1790 
1800 
1810 
1820 
1830 
1840 
1860 
1860 
1870 
2000 
2010 
2020 
2030 
2040 
2060 
2060 
2070 
2080 


FOR  I!-l  TO  3:P0S2(I!)-P0S1(I!)*CD+CVEC1(I!)*SD:NEXT 
L2=»FNATN2(P0S2(2),P0S2(D) :P2-FNASN(P0S2(3)) : P-P2 : L-L2 : G08UB  1760 
X-GVEC(1)*EA8T(1)+GVEC(2)*EA8T(2)+GVEC(3)*EAST(3) 
Y=-(GVEC(1)*N0RTH(1)+GVEC(2)*N0RTH(2)+GVEC(3)*N0RTH(3)) 

A21-FNATNP(Y,X):L2— L2: RETURN :* CONVERT  BACK  TO  LEFT-HANDED  OUTPUT 

• 

'INVER8E  SOLUTION,  SPHER  EARTH- RECT  COORD.  ALL  ANGLES  MUST  BE  IN  RADIANS 

•INPUT:  LATITUDES  PI  ft  P2,  AND  LONGITUDES  LI  ft  L2. 

•OUTPUT:  DISTANCE  DD  TO  A  POINT  P2.   (NOTE:  0  <=  DD  <=  PI  RADIAN8) . 

'  FORWARD  AZIMUTH  A12.  AND  BACKWARD  AZIMUTH  A21. 

P=P1:L=-L1: GOSUB  1750:'  CHNG  SIGN  OF  LI  GIVES  RIGHT-HANDED  COORDS 

FOR  I!-l  TO  3:P0S1(I!)-P0SI(I!):N0RTH1(I!)-N0RTH(I!):EA8T1(I!)-EAST(I!): 

NEXT 

P=P2:L=-L2: GOSUB  1750:'     CHNG  SIGN  OF  L2  GIVES  RIGHT-HANDED  COORDS 

FOR  I!-l  TO  3:P0S2(I!)-P0SI(I!):N0RTH2(I!)-N0RTH(I!):EAST2(I!)-EAST(I!): 

NEXT 

DD-FNACS(P081(1)*P082(1)+P0S1(2)*P0S2(2)+P0S1(3)*P0S2(3)) 

X=FNCX(P0glU).P0Sl(2)JP0Sl(3).P0S2(l),P0S2(2)1P0S2(3)) 

Y-FNCY(POSKl)  ,P0S1(2)  ,P0S1  (3)  ,P0S2(1)  ,P0S2(2)  ,P0S2(3)) 

Z-FNCZ(P081(1),P081(2)>P0S1(3)IP0S2(1),P0S2(2)JP0S2(3)) 

A12-FNATNP(X*NORTH1(1)+Y*NORTH1(2)+Z*NORTH1(3)J 

-(X*EAST1(1)*Y*EAST1(2)+Z*EAST1(3))) 

A21-FNATNP(-(X*N0RTH2(1)+Y*N0RTH2(2)+Z*N0RTH2(3)), 

X*EA8T2(1)+Y*EAST2(2)+Z*EAST2(3)) 

RETURN 

• 

CA-COS (A) : SA-SIN(A) : >Y*CA+Z*SA : Z«Z*CA-Y*SA : Y-T : RETURN : '  X-AXIS  ROT 
CA=COS (A) :SA=SIN(A):T=Z*CA+X*SA:X=X*CA-Z*SA:Z=T: RETURN:'  Y-AXIS  ROT 

CA»C0S(A):SA-SIN(A):T-X*CA+Y*8A:Y-Y*CA-X*SA:X-T:RETURN:*  Z-AXIS  ROT 

• 

'UNIT  VECTORS:  POSITION,  NORTH  ft  EAST. 

SL=SIN(L) :CL=C08(L) :SP=SIN(P) :CP=COS(P) 

P0SI(1)=CP*CL:P08I(2)=<:P*SL:P0SI(3)=SP 

N0RTH(1)=-SP*CL: NORTH (2)=-SP*SL:N0RTH(3)=CP 

EASTC  D—SL :  EAST(2)-CL :  EAST(3)-0 :  RETURN 
• 

' VECTORS :CVEC=COURSE  ft  GVEC=GREAT  CIRCLE  NORMAL 

X=0:Y=0:Z=1:A=A12:GOSUB  1700:A=P:G08UB  1710:A=-L:G0SUB  1720 

CVEC1(1)-X:CVEC1(2)-Y:CVEC1(3)-Z 

GVEC(1)-FNCX(P0S1(1),P081(2),P0S1(3),CVEC1(1),CVEC1(2),CVEC1(3)) 

GVEC(2)»FNCY(P081(1),P0S1(2),P081(3)<CVEC1(1),CVEC1(2),CVEC1(3)) 

GVEC(3)aFNCZ(P0Sl(l),P0Sl(2),P081(3),CVECl(l),CVECl(2),CVECl(3)) 

RETURN 

CLS: PRINT  SPC(20); "ALGORITHM  DEMO 

PRINT: PRINT: PRINT 

PRINT  SPC(15);"1)  DIRECT  SOLUTION 

-2)  INVERSE  SOLUTION 

"3)  FIND  CPA 

■4)  SPEED  NEEDED  TO  INTERCEPT 

"5)  TIME  NEEDED  TO  INTERCEPT 
PRINT:PRINT  SPC(15);"6)  qUIT 
GOSUB  9010:C-VAL(C$):ON  C  GOSUB  3000,4000,5000,6000,7000,8000 


PRINT  SPC(15) 
PRINT  SPC(15) 
PRINT  SPCC15) 
PRINT  SPC(15) 
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2090  GOTO  2000 

2100  ' 

3000  CLS:PRINT  8PC(15); "DIRECT  SOLUTION" : PRINT: PRINT 


PRINT  "1st  LATITUDE  DD.MM8S  (-8)  " 

1150:P1-X*RD 

PRINT  "lat  LONGITUDE  DDD.MMSS  (-E) 

1150:L1-X*RD 

PRINT  "INITIAL  COURSE  DDD.MMSS  "; 

1150:A12-X*RD 

INPUT  "DISTANCE  (n.  Mi.)  ?  ",D 


3010  PRINT  SPC(IO); 

3020  INPUT  V$:G08UB 

3030  PRINT  SPC(IO); 

3040  INPUT  V|: GOSUB 

3060  PRINT  SPC(IO); 

3060  INPUT  V$: GOSUB 

3070  PRINT  SPC(IO); 

3080  Dl-D*RD/80# 

3000  MO-100 

3100  DD-D1:G0SUB  1340: PRINT: PRINT  SPC (8) ; "SPHERICAL  EARTH  DIRECT  SOLUTION 

3110  PRINT  SPC(12);"2nd  LATITUDE  "; :X«P2/RD: GOSUB  1090:PRINT  V$ 

3120  PRINT  SPC(12);"2nd  LONGITUDE  " ; : X-L2/RD : GOSUB  1000: PRINT  V$ 

3130  PRINT  SPC(12) ; "BACK  AZIMUTH  ■; :X«A21/RD: GOSUB  1000:PRINT  V$ 

3140  GOSUB  9000: GOTO  2000 

3150  * 

4000  CLS: PRINT  SPC(15); "INVERSE  SOLUTION": PRINT: PRINT 


4010  PRINT  SPC(IO); 

4020  INPUT  VI: GOSUB 

4030  PRINT  SPC(IO); 

4040  INPUT  V$: GOSUB 

4060  PRINT  SPC (10); 

4060  INPUT  V$: GOSUB 

4070  PRINT  SPC(IO); 


PRINT  "1st  LATITUDE  DD.MMS8  (-8)  "; 

1160.P1-X 

PRINT  "1st  LONGITUDE  DDD.MM88  (-E)  "; 

1150:L1-X 

PRINT  "2nd  LATITUDE  DD.MMSS  (-8)  •; 

1160:P2-X 

: PRINT  "2nd  LONGITUDE  DDD.MMSS  (-E)  "; 
4080  INPUT  VI  GOSUB  1150:L2-X 
4090  P1-P1+RD : P2»P2*RD : L1-L1*RD : L2-L2*RD 
4100  DD-D1: GOSUB  1640 

4110  PRINT:PRINT  SPC (8) ; "SPHERICAL  INVERSE  SOLUTION 
4120  ' 

4130  PRINT  SPC(12); "DISTANCE  "; 
4140  MO100:  PRINT  FNR(60#*DD/RD) ; ■  n.ml. 

4150  PRINT  SPC(12) ; "FORWARD  COURSE  ■; :X-A12/RD: GOSUB  1090:PRINT  V$ 
4160  PRINT  SPC ( 12 ); "BACK  COURSE  ■; :X»A21/RD: GOSUB  1090:PRINT  V$ 
4170  GOSUB  9000: GOTO  2000 
4180  ' 

5000  CLS: PRINT  SPC(20);"FIND  CPA" : PRINT : PRINT 
5010  PRINT  8PC( 10);: PRINT  "lat  LATITUDE  DD.MM8S  (-8)  ■; 
5020  INPUT  V$: GOSUB  1150:P1-X*RD 

5030  PRINT  SPC (10);: PRINT  "lat  LONGITUDE  DDD.MM8S  (-E)  "; 
5040  INPUT  V$: GOSUB  1150:L1-X*RD 

5060  PRINT  8PC(10);: PRINT  "INITIAL  COURSE  DDD.MMSS  "; 
5060  INPUT  V): GOSUB  1150:A1»X*RD 
5070  PRINT  8PC(10);: INPUT  "SPEED  (knota)  ?  ",S1 
5080  PRINT  8PC( 10);: PRINT  "2nd  LATITUDE  DD.MM8S  (-8)  "; 
5090  INPUT  V$: GOSUB  1160:P2=X*RD 

5100  PRINT  8PC( 10);: PRINT  "2nd  LONGITUDE  DDD.MMSS  (-E)  •; 
5110  INPUT  V$ : GOSUB  11S0:L2-X*RD 

5120  PRINT  8PC( 10);: PRINT  "INITIAL  COURSE  DDD.MMSS  "; 
5130  INPUT  V$: GOSUB  1160:A2-X*RD 
5140  PRINT  SPC(IO);: INPUT  "SPEED  (knota)  ?  ",82 
5160  Bl-Sl/AE:B2-82/AE 
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5160  ' 

5170  P=P1:L=-L1:G0SUB  1750:'  CHNG  SIGN  OF  LI  GIVES  RIGHT-HANDED  C00RD8 

5180  FOR  Il-l  TO  3:X1(I!)-P08I(I!):NEXT 

5190  X=0:Y=0:Z=1:A=A1:GOSUB  1700 :A=P: GOSUB  1710:A=-L:G08UB  1720 

5200  C1(1)"X'C1(2)"Y'C1(3)"Z 

5210  P=P2;L=-L2: GOSUB  1750:'   CHNG  SIGN  OF  L2  GIVES  RIGHT-HANDED  COORDS 

5220  FOR  Il-l  TO  3:X2(I!)-P0SI(I!) :NEXT 

5230  X»0:Y-0:Z«1:A-A2: GOSUB  1700 :A»P: GOSUB  1710: A— L:GOSUB  1720 

5240  C2(1)-X:C2(2)-Y:C2(3)-Z 

5260  ' 

5260  X1C2-X1(1)*C2(1)+X1(2)*C2(2)+X1(3)*C2(3) 

5270  C1X2=C1(1)*X2(1)+C1(2)*X2(2)+C1(3)*X2(3) 

5280  X1X2-X1(1)*X2(1)+X1(2)*X2(2)+X1(3)*X2(3) 

5200  C1C2=C1(1)*C2(1)+C1(2)*C2(2)+C1(3)*C2(3) 

5300  BA=-B1*X1C2  -  B2*C1X2 

5310  BB-  B1*C1X2  t  B2*X1C2 

6320  BC»-B1*X1X2  ♦  B2*C1C2 

5330  BD-  B1*C1C2  -  B2*X1X2 

5340  ' 

5360  T=l : IT ! =0 : '  ITERATE  WITH  NEWTON-RAPHSON 

5360  B1T-B1*T:B2T-B2*T:S1-SIN(B1T)  :C1-C0S(BIT)  S2-SIN(B2T)  :C2-C0S(B2T) 

5370  S1S2-S1*S2:C1C2-C1*C2;S1C2»S1*C2:C1S2=C1*S2 

5380  F=BA*S1S2+BB*C1C2+BC*S1C2+BD*C1S2 

5300  FP-- (BC*B2+BD*B1)*S182+(BD*B2+BC*B1)*C1C2+(BA*B2-BB*B1)*S1C2- 

(BB*B2-BA*B1)*C1S2 

5400  IT!-IT1+1:C0RR-F/FP:T-T-C0RR:IF  ABS(CORR)<  00001  THEN  5440 

5410  IF  IT!>50  THEN  PRINT  "NO  CONVERGENCE* : STOP 

6420  GOTO  6360 

5430  ' 

5440  B1T-B1*T:B2T-B2*T:CB1T»C0S(B1T):SB1>8IN(B1T):CB2T-C0S(B2T):SB2T-SIN(B2T) 

6460  FOR  I!-l  TO  3:P0S1(I!)-X1(I!)*CB1T+C1(I!)*SB1T 

5460  P0S2(I!)-X2(I!)*CB2TtC2(I!)*SB2T:NEXT  I! 

5470  P1=FNASN(P0S1(3)):L1=-FNATN2(P081(2),P0S1(1)) 

5480  P2-FNASN(PQS2(3)):L2— FNATN2(P0S2(2)(P0S2(1)) 

5400  GO SUB  1540 

5500  X-T:GOSUB  1010: PRINT: PRINT  SPC(10);"TINE  TO  CPA  ■  ";V$ 

5510  M0»100#:PRINT  SPC(IO); "DISTANCE  AT  CPA  ■  ";FNR(60#*DD/RD);"  n.mi. 

5620  PRINT  SPC(IO); "BEARING  AT  CPA  -  ■ ; :X»A12/RD:G0SUB  1000: PRINT  V$ 

5530  PRINT  SPC(10);"NO.  ITERATIONS  -  ";IT! 

5540  TCPA=T:GOSUB  9000 

5660  ' 

5560  CLS: PRINT  SPC(22);"FIND  CPA": PRINT 

5670  PRINT  "  TIME  LAT  1  LONG  1  DISTANCE  BEARING(l->2)" 

5580  D>TCPA/6#:T-0:FOR  Tl-l  TO  13:B1T-B1*T:B2T-B2*T 

5590  FOR  I!-l  TO  3:P081(I!)-X1(I!)*C0S(B1T)+C1(I!)*SIN(B1T) 

5600  P0S2(I!)-X2(I!)*C0S(B2T)+C2(I!)*SIN(B2T):NEXT  II 

5610  P1=FNASN(P081(3)):L1=-FNATN2(P0S1(2).P081(1)) 

5620  P2-FNASN(P0S2(3)):L2— FNATN2(P0S2(2),P0S2(1)) 

6630  GOSUB  1540 

5640  X-T: GOSUB  1010 .PRINT  V$; 

5660  LOCATE  CSRLIN, 12: X-P1/RD: GOSUB  1090: PRINT  V$; 

5660  LOCATE  CSRLIN, 23 :X-L1/RD: GOSUB  1090: PRINT  V$; 

6670  LOCATE  CSRLIN , 37 : PRINT  FNR(60#*DD/RD) ; 
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5680  LOCATE  CSRLIN,52:X-A12/RD  GOSUB  1000: PRINT  V$ 

6690  T-T+DT:NEXT  T! 

5700  GOSUB  9000: GOTO  2000 

5710  ■ 

6000  CL8: PRINT  SPC ( 10) ;" SPEED  NEEDED  TO  INTERCEPT  (l->2)" : PRINT: PRINT 

6010  CLNX-O: GOSUB  6610:'  GET  INPUT 

6020  CLNX-10: LOCATE  CLNX,1:FQR  H-l  TO  11: PRINT  8PC(79) :NEXT: '  CLEAR  SCREEN 

6030  LOCATE  CLN% , 1 1 : PRINT  "TIME  TO  INTERCEPT  (SECONDS)  ";: INPUT  TMI 

6040  TMI-TMI/3600# : TM-TMI : PI-PIS : L1-L1S : P2-P2S : L2-L2S : A2-A28 

6060  ' 

6060  '  COMPUTE  SPEED 

6070  P=P18:L=-L1S: GOSUB  1750:'  CHNG  SIGN  OF  LI  GIVES  RIGHT-HANDED  C00RD8 

6080  FOR  I!-l  TO  3:X1(I!)«P08I(I!) :NEXT 

6090  P»P2S:L— L2S: GOSUB  1750:'  CHNG  SIGN  OF  L2  GIVES  RIGHT-HANDED  COORDS 

6100  FOR  Ii-1  TO  3:X2(I!)-P0SI(I!):NEXT 

6110  X»0:Y»0:Z=1:A*A2: GOSUB  1700  :A=P:  GOSUB  1710  :A—L:  GOSUB  1720 

6120  C2(1)-X:C2(2)=Y:C2(3)»Z 

6130  B2T=B2*TM:CB2T=C0S(B2T) :SB2T=SIN(B2T) :CS=0# 

6140  FOR  Il-l  TO  3:P082(I!)-X2(I!)*CB2T+C2(I!)*SB2T 

6160  CS«CS+X1(I!)*P0S2(I!):NEXT  I! 

6160  S-FNACS(CS) :SPD-8*AE/(TM-EPS*(TM-0)) 

6170  " 

6180  '  GET  INVERSE  SOLN 

6190  P2«FNASN(P0S2(3)):L2— FNATN2(P082(2) ,P0S2(1)):G08UB  1640 

6200  ' 

6210  MO-10: LOCATE  CLNX+2, 11: PRINT  "SPEED  REQUIRED  -  ";FNR(SPD);"  knots" 

6220  GOSUB  6810: '  PRINT  OUT  BEARING,  RANGE,  LAT  k  LONG 

6230  LOCATE  CLNX+8, 15: PRINT  "1)  CHANGE  TIME  OF  INTERCEPT" 

6240  LOCATE  CLNX+O, 15: PRINT  "2)  NEW  PROBLEM" 

6250  LOCATE  CLNt+10, 16:PRINT  "3)  MASTER  MENU" 

6260  GOSUB  9010:C»VAL(C$) :0N  C  GOTO  6020,6000,2000 

6270  GOTO  6260 

6280  ' 

8500  '  INPUT  ROUTINE 

6610  LOCATE  CLNX+3 , 1 1 : PRINT  "let  LATITUDE  DD.MMSS  (-S)  ■; 

6620  INPUT  V$: GOSUB  1150:P1S=X*RD 

6530  LOCATE  CLNX+4 , 1 1 : PRINT  "1st  LONGITUDE  DDD.MMSS  (-E)  "; 

6640  INPUT  V$: GOSUB  1160:L1S-X*RD 

6560  LOCATE  CLNX+6, 11: PRINT  "2nd  LATITUDE  DD.MM8S  (-S)  ■; 

6660  INPUT  V$: GOSUB  1160:P2S=X*RD 

6570  LOCATE  CLNX+6, 11: PRINT  "2nd  LONGITUDE  DDD.MMSS  (-E)  ■; 

6580  INPUT  V$: GOSUB  1160:L2S~X*RD 

6590  LOCATE  CLNX+7, 11: PRINT  "2nd  COURSE  DDD.MMSS  ■; 

6600  INPUT  V$: GOSUB  1150:A2S-X*RD 

6610  LOCATE  CLNX+8, 11: INPUT  "2nd  SPEED  (knots)  ?  ",S2 

6620  B2-S2/AE: RETURN 

6630  ' 

6800  '  OUTPUT  BEARING,  RANGE,  LAT  *  LONG 

6810  LOCATE  CLNX+3 , 1 1 : PRINT  "BEARING  TO  INTERCEPT  =  "; 

6820  X=A12/RD:G08UB  1090: PRINT  V$ 

6830  LOCATE  CLNX+4 . 1 1 : PRINT  "RANGE  TO  INTERCEPT  =  ";FNR(60#*DD/RD) ; "  n.ai." 

6840  LOCATE  CLNX+5, 11: PRINT  "LATITUDE  OF  INTERCEPT  =  "; 

6860  X=P2/RD: GOSUB  1090: PRINT  V$ 
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8860  LOCATE  CLNl+6 , 1 1 : PRINT  "LONGITUDE  OF  INTERCEPT  -  "; 

8870  X-L2/RD: GOSUB  1000: PRINT  V$ 

8880  RETURN 

6890  ' 

7000  CLS:PRINT  SPC(10);"TIME  NEEDED  TO  INTERCEPT  (l->2)": PRINT: PRINT 

7010  CLNX-O: GOSUB  6610:'  GET  INPUT 

7020  • 

7030  *  FIND  Vain  AND  Tain  vel. 

7040  '  Compute  arc  distance  S 

7060  P=P1S:L=-L1S:G0SUB  1750:*  CHNG  SIGN  OF  LI  GIVES  RIGHT-HANDED  COORDS 

7060  FOR  I!=l  TO  3:11(1! )=P08I(I!)  NEXT 

7070  P=P2S;L=-L2S:G0SUB  1750:'  CHNG  SIGN  OF  L2  GIVES  RIGHT-HANDED  COORDS 

7080  FOR  IJ-1  TO  3:X2(I!)=P08I(I!) :NEXT 

7090  X=0:Y=0:Z=1:A=A2S: GOSUB  1700 :A=P: GOSUB  1710 :A=-L: GOSUB  1720 

7100  C2(1)-X:C2(2)-Y:C2(3)-Z 

7110  ' 

7120  ■  INITIALIZE 

7130  X1X2-X1(1)*X2(1)+X1(2)*X2(2)+X1(3)*X2(3) 

7140  X1C2-X1(1)*C2(1)+X1(2)*C2(2)+X1(3)*C2(3) 

7150  • 

7160  '  BEGIN  ITERATION 

7170  S-FNACS (X1X2) : T-S+AE/S2 : IT ! "1 

7180  B2T-B2*T:CB2T»CQS(B2T)  :SB2T-SIN(B2T) 

7190  CS-X1X2*CB2T+X1C2*SB2T: S-FNACS (CS) : SS-SIN(S) 

7200  DSDT=(X1X2*SB2T-X1C2*CB2T) *B2/SS 

7210  F=T*DSDT-8 

7220  FP-T*(B2*B2*(X1X2*CB2T+X1C2*SB2T)-CS*DSDT*DSDT)/SS 

7230  IT!sIT!+l:CORR=F/FP:T=T-CORR:IF  ABS(CORR)<. 00001  THEN  7270 

7240  IF  IT!>50  THEN  PRINT  "NO  CONVERGENCE" : STOP 

7250  GOTO  7180 

7260  ' 

7270  TM8-T:VMIN-S*AE/T 

7280  LOCATE  CLNX+10,11:M0-10 

7290  PRINT  "MINIMUM  SPEED  REQUIRED  TO  INTERCEPT  -  ";FNR(VMIN) ;"  knots" 

7300  LOCATE  CLNX+11,11:X-TMS:G0SUB  1010 

7310  PRINT  "TIME  REQ/D  TO  INTERCEPT  AT  MIN  SPEED  -  ";V$ 

7320  ' 

7330  CLNX-13: LOCATE  CLNt,l:FOR  Il-l  TO  11:PRINT  SPC(79) :NEXT: *  CLEAR  SCREEN 

7340  LOCATE  CLN%,11:PRINT  "INTERCEPTOR  SPEED  (knots)  ";:INPUT  SPD 

7360  IF  8PD>=VMIN  THEN  7390 

7360  LOCATE  CLNX+2 , 1 1 : PRINT  "SPEED  TOO  LOW,  CANNOT  INTERCEPT" 

7370  GOSUB  9000: GOTO  7330 

7380  ' 

7390  Bl-SPD/AE:T-TM8/5:  T-.l  :IT!-1 

7400  B2T»B2*T:CB2T-C0S(B2T) :SB2T»8IN(B2T) 

7410  CS-X1X2*CB2T+X1C2*SB2T: S-FNACS (CS) : SS-SIN(S) 

7420  DSDT=(X1X2*SB2T-X1C2*CB2T) *B2/SS 

7430  F-S/T-B1:FP-(DSDT-S/T)/T 

7440  IT!=IT!+1:C0RR*F/FP:T»T-C0RR:IF  ABS(CORR)<. 00001  THEN  7500 

7460  IF  IT!>60  THEN  PRINT  "NO  CONVERGENCE" : STOP 

7460  IF  ABS(CORR)< 1000000000*  THEN  7400 

7470  LOCATE  CLNX+2, 11: PRINT  "SPEED  TOO  HIGH,  NO  CONVERGENCE" 

7480  LOCATE  CLNX+4: GOSUB  9000: GOTO  7330 
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7400  ' 

7500  X=T:G08UB  1010 

7510  LOCATE  CLN%+2 , 1 1 : PRINT  "TIME  REQUIRED  -  ";V$ 

7520  ' 

7530  '  GET  INVERSE  SOLN 

7540  FOR  Il-l  TO  3:P0S2(I!)-X2(I!)*CB2T+C2(I!)*SB2T:NEXT  I! 

7550  P1=P18:L1=L1S:P2=FNASN(P0S2(3)):L2=-FNATN2(P0S2(2)1P0S2(1)) :GOSUB  1540 

7560  ■ 

7570  GOSUB  6810: '  PRINT  OUT  BEARING,  RANGE,  LAT  *  LONG 

7580  LOCATE  CLNX+8, 15: PRINT  "1)  CHANGE  INTERCEPTOR  SPEED" 

7500  LOCATE  CLNX+O, 15: PRINT  "2)  NEW  PROBLEM" 

7600  LOCATE  CLNt+10, 15: PRINT  "3)  MASTER  MENU" 

7610  GOSUB  0010:C-VAL(C$):ON  C  GOTO  7330,7000,2000 

7620  GOTO  7610 

7630  ' 

8000  CL3:END 

8010  * 

0000  PRINT: PRINT  SPC(IO); "PRESS  ANY  KEY  TO  CONTINUE 

0010  FOR  I!=l  TO  0:C$=INKEY$:NEXT 

0020  C$=INKEY$:IF  C$-"  THEN  0020 

0030  RETURN 
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